package scu.maqiang.mesh;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;
import java.util.*;
import java.util.stream.IntStream;

//import com.sun.tools.javac.code.Attribute.Array;

import scu.maqiang.fes.EleFunc;
import scu.maqiang.numeric.*;

import static java.lang.Math.*;


/**
 * 
 * @author Administrator
 *
 */
public class Mesh2T3 extends Mesh {

	public Mesh2T3() {
		super();
		nDim = 2;
		nPerElement = 3;
		nPerBoundary = 2;
		nBoundaryPerElement = 3;
		nEdgePerElement = 3;
		nEdgePerBoundary = 1;
		tecplotType = "FETRIANGLE";
		meditMeshType = new String[]{"Vertices", "Edges", "Points"};
		gmshType = new String[]{"ST", "VT", "TT"};
		gmshElementType = 2;
	}

	@Override
	public boolean checkMesh(double[][] coord) {
		return CG.area2DT3(coord) > 0;
	}

	public int[] pointInMesh(double[] xy) {
		this.sortElements();
		int left = 0;
		int right = nt - 1;
		int middle = 0;
		double[][] coord = new double[nPerElement][nDim];
		double[] N = new double[nPerElement];
		int[] ele = null;
		int size = 5;
		while(left <= right) {
			middle = (left + right) / 2;
			ele = elements[middle];
			getCoordPerBoundary(ele, coord);
			double[] p0 = getinteriorPointInElement(coord);
			if (xy[0] < p0[0]) {
				right = middle;
			} else {
				left = middle;
			}
			if(right - left < size) {
				break;
			}
		}
		
		System.out.println("left:" + left + "\tright: " + right);
		for(int i = Math.max(0, left - 10); i <= Math.min(right + 10, nt); i++) {
			ele = elements[i];
			getCoordPerElement(ele, coord);
			double[] p0 = getinteriorPointInElement(coord);
			System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
			//System.out.println(MVO.toString(coord));
			//CG.T3ShapeFunction(coord, xy, N);
			//System.out.println(MVO.toString(N));
			if(CG.pointInT3(coord, xy, N)) {
				//System.out.println(Arrays.toString(N));
				//System.out.println(Arrays.toString(ele) + "\t" + Arrays.toString(p0));
				return ele;
			}
		}
		return null;
	}

	public int[] signsInEdgesByVertex(int[] ele) {
		int[] signs = new int[3];
		signs[0] = ele[0] < ele[1] ? 1:-1;
		signs[1] = ele[1] < ele[2] ? 1:-1;
		signs[2] = ele[2] < ele[0] ? 1:-1;
		return signs;
	}

	public int[][] signsInEdgesByElement(SIMatrix edgeMat) {
		int[][] signs = new int[nt][3];
		for(int i = 0; i < nt; i++) {
			int[] ele = getElement(i);
			int ii = edgeMat.getElement(ele[1], ele[0]);
			signs[i][0] = ii == 0?1:(i < ii?1:-1);
			ii = edgeMat.getElement(ele[2], ele[1]);
			signs[i][1] = ii == 0?1:(i < ii?1:-1);
			ii = edgeMat.getElement(ele[0], ele[2]);
			signs[i][2] = ii == 0?1:(i < ii?1:-1);
		}
		return signs;
	}

	public Mesh2T3 refine() {
		Mesh2T3 newMesh = new Mesh2T3();
		SIMatrix C = this.renumberingEdges();
		System.out.println(C);
		int nEdge = C.getNNZ();
		int newNv = this.nv + nEdge;
		int newNt = this.nt * 4;
		int newNb = this.nb * 2;
		newMesh.initMesh0(newNv, newNt, newNb);
		for(int i = 0; i < nv; i++) {
			newMesh.nodes[i] = this.nodes[i].clone();
			newMesh.nodeLabel[i] = this.nodeLabel[i];
		}
		int[] edgI = new int[nEdge];
		int[] edgJ = new int[nEdge];
		int[] edgeIdx = new int[nEdge];
		C.getTriValues(edgI, edgJ, edgeIdx);
		for(int i = 0; i < nEdge; i++) {
			double[] node1 = this.getNode(edgI[i]);
			double[] node2 = this.getNode(edgJ[i]);
			newMesh.nodes[nv + i] = new double[] {0.5 * (node1[0] + node2[0]), 0.5 * (node1[1] + node2[1])};
		}
		//newMesh.displayNodes();
		int eleIdx = 0;
		for(int i = 0; i < nt; i++) {
			int[] ele = elements[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + nv - 1;
			int jk = Math.max(C.getElement(ele[1], ele[2]), C.getElement(ele[2], ele[1])) + nv - 1;
			int ki = Math.max(C.getElement(ele[2], ele[0]), C.getElement(ele[0], ele[2])) + nv - 1;
			newMesh.elements[eleIdx] = new int[]{ele[0], ij, ki};
			newMesh.elementLabel[eleIdx++] = this.elementLabel[i];
			newMesh.elements[eleIdx] = new int[]{ele[1], jk, ij};
			newMesh.elementLabel[eleIdx++] = this.elementLabel[i];
			newMesh.elements[eleIdx] = new int[]{ele[2], ki, jk};
			newMesh.elementLabel[eleIdx++] = this.elementLabel[i];
			newMesh.elements[eleIdx] = new int[]{ij, jk, ki};
			newMesh.elementLabel[eleIdx++] = this.elementLabel[i];
		}
		//newMesh.displayElements();
		eleIdx = 0;
		for(int i = 0; i < nb; i++) {
			int[] ele = boundaries[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + nv - 1;
			newMesh.boundaries[eleIdx] = new int[]{ele[0], ij};
			newMesh.boundaryLabel[eleIdx++] = this.boundaryLabel[i];
			newMesh.boundaries[eleIdx] = new int[]{ij, ele[1]};
			newMesh.boundaryLabel[eleIdx++] = this.boundaryLabel[i];
		}
		newMesh.displayBoundarys();
		return newMesh;
	}

	public Mesh2T3 selfRefine() {
		SIMatrix C = this.renumberingEdges();
		System.out.println(C);
		int nEdge = C.getNNZ();
		int oldNv = nv;
		int oldNt = nt;
		int oldNb = nb;
		int newNv = this.nv + nEdge;
		int newNt = this.nt * 4;
		int newNb = this.nb * 2;
		double[][] oldNodes = this.nodes;
		int[] oldeNodeLabel = this.nodeLabel;
		int[][] oldElements = this.elements;
		int[] oldElementLabel = this.elementLabel;
		int[][] oldBoundaries = this.boundaries;
		int[] oldBounadryLabel = this.boundaryLabel;
		initMesh0(newNv, newNt, newNb);
		//System.out.println(newNv + MVO.toString(oldNodes));
		for(int i = 0; i < oldNv; i++) {
			nodes[i] = oldNodes[i];
			nodeLabel[i] = oldeNodeLabel[i];
		}
		int[] edgI = new int[nEdge];
		int[] edgJ = new int[nEdge];
		int[] edgeIdx = new int[nEdge];
		C.getTriValues(edgI, edgJ, edgeIdx);
		for(int i = 0; i < nEdge; i++) {
			double[] node1 = nodes[edgI[i]];
			double[] node2 = nodes[edgJ[i]];
			nodes[oldNv + i] = new double[] {0.5 * (node1[0] + node2[0]), 0.5 * (node1[1] + node2[1])};
		}
		displayNodes();
		int eleIdx = 0;
		System.out.println(oldNt);
		for(int i = 0; i < oldNt; i++) {
			int[] ele = oldElements[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + oldNv - 1;
			int jk = Math.max(C.getElement(ele[1], ele[2]), C.getElement(ele[2], ele[1])) + oldNv - 1;
			int ki = Math.max(C.getElement(ele[2], ele[0]), C.getElement(ele[0], ele[2])) + oldNv - 1;
			elements[eleIdx] = new int[]{ele[0], ij, ki};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ele[1], jk, ij};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ele[2], ki, jk};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ij, jk, ki};
			elementLabel[eleIdx++] = oldElementLabel[i];
		}
		displayElements();
		eleIdx = 0;
		for(int i = 0; i < oldNb; i++) {
			int[] ele = oldBoundaries[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + oldNv - 1;
			boundaries[eleIdx] = new int[]{ele[0], ij};
			boundaryLabel[eleIdx++] = oldBounadryLabel[i];
			boundaries[eleIdx] = new int[]{ij, ele[1]};
			boundaryLabel[eleIdx++] = oldBounadryLabel[i];
		}
		return this;
	}

	public Mesh2T3 selfRefine(VectorFunc bdFunc, double[][] param) {
		SIMatrix C = this.renumberingEdges();
		System.out.println(C);
		int nEdge = C.getNNZ();
		int oldNv = nv;
		int oldNt = nt;
		int oldNb = nb;
		int newNv = this.nv + nEdge;
		int newNt = this.nt * 4;
		int newNb = this.nb * 2;
		double[][] oldNodes = this.nodes;
		int[] oldeNodeLabel = this.nodeLabel;
		int[][] oldElements = this.elements;
		int[] oldElementLabel = this.elementLabel;
		int[][] oldBoundaries = this.boundaries;
		int[] oldBounadryLabel = this.boundaryLabel;
		initMesh0(newNv, newNt, newNb);
		//System.out.println(newNv + MVO.toString(oldNodes));
		for(int i = 0; i < oldNv; i++) {
			nodes[i] = oldNodes[i];
			nodeLabel[i] = oldeNodeLabel[i];
		}
		int[] edgI = new int[nEdge];
		int[] edgJ = new int[nEdge];
		int[] edgeIdx = new int[nEdge];
		C.getTriValues(edgI, edgJ, edgeIdx);
		for(int i = 0; i < nEdge; i++) {
			double[] node1 = nodes[edgI[i]];
			double[] node2 = nodes[edgJ[i]];
			nodes[oldNv + i] = new double[] {0.5 * (node1[0] + node2[0]), 0.5 * (node1[1] + node2[1])};
		}
		displayNodes();
		int eleIdx = 0;
		System.out.println(oldNt);
		for(int i = 0; i < oldNt; i++) {
			int[] ele = oldElements[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + oldNv - 1;
			int jk = Math.max(C.getElement(ele[1], ele[2]), C.getElement(ele[2], ele[1])) + oldNv - 1;
			int ki = Math.max(C.getElement(ele[2], ele[0]), C.getElement(ele[0], ele[2])) + oldNv - 1;
			elements[eleIdx] = new int[]{ele[0], ij, ki};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ele[1], jk, ij};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ele[2], ki, jk};
			elementLabel[eleIdx++] = oldElementLabel[i];
			elements[eleIdx] = new int[]{ij, jk, ki};
			elementLabel[eleIdx++] = oldElementLabel[i];
		}
		displayElements();
		eleIdx = 0;
		for(int i = 0; i < oldNb; i++) {
			int[] ele = oldBoundaries[i];
			int ij = Math.max(C.getElement(ele[0], ele[1]), C.getElement(ele[1], ele[0])) + oldNv - 1;
			if (bdFunc != null) {
				nodes[ij] = bdFunc.action(getNode(ij), oldBounadryLabel[i], param);
			}
			boundaries[eleIdx] = new int[]{ele[0], ij};
			boundaryLabel[eleIdx++] = oldBounadryLabel[i];
			boundaries[eleIdx] = new int[]{ij, ele[1]};
			boundaryLabel[eleIdx++] = oldBounadryLabel[i];
		}
		return this;
	}

	public void selfRefine(int[] elements, ScalarFunc bdFunc, double[] param) {
		int[] refineIdx = new int[nv];
		IntStream.of(elements).forEach(ele -> refineIdx[ele] = 3);
	}
	
	/**
	 * ��������ӵ�ԭ������
	 *
	 *            ������
	 * @return ԭ����������ʵ������ĳ����ϲ�
	 */
//	public Mesh2DT3 add(Mesh2DT3 meshB) {
//		int[] newNodesB, newBdB, ABd;
//		double tol = 1.0e-3, xA, yA, xB, yB;
//		int i, j, nvA, nvB, nConnected, n, nBd;
//		int[] idx = new int[2];
//		int[] idx2 = new int[2];
//		int[] temp = new int[2];
//		newNodesB = new int[meshB.nv];
//		newBdB = new int[meshB.nb];
//		ABd = new int[nb];
//		Arrays.fill(newNodesB, -1);
//		Arrays.fill(newBdB, -1);
//		Arrays.fill(ABd, -1);
//		nvA = nv;
//		nvB = meshB.nv;
//		for (i = 0; i < nvA; i++) {
//			xA = nodes[i][0];
//			yA = nodes[i][1];
//			for (j = 0; j < nvB; j++) {
//				xB = meshB.nodes.get(j).coord[0];
//				yB = meshB.nodes.get(j).coord[1];
//				if (abs(xA - xB) < tol && abs(yA - yB) < tol) {
//					newNodesB[j] = i;
//					break;
//				}
//			}
//		}
//
//		for (i = 0; i < meshB.nb; i++) {
//			idx = meshB.boundarys.get(i).idx;
//			temp[0] = newNodesB[idx[0]];
//			temp[1] = newNodesB[idx[1]];
//			if (temp[0] > -1 && temp[1] > -1) {
//				newBdB[i] = 1;
//				for (j = 1; j < nb; j++) {
//					idx2 = boundarys.get(j).idx;
//					if ((idx2[0] == temp[0] && idx2[1] == temp[1]) || (idx2[1] == temp[0] && idx2[1] == temp[0])) {
//						ABd[j] = 1;
//						break;
//					}
//				}
//			}
//		}
//
//		// nConnected = (int)(IntStream.of(newNodesB).filter(iidx -> iidx !=
//		// -1).count());
//		// nBd = (int)(IntStream.of(newBdB).filter(iidx -> iidx != -1).count());
//		nConnected = 0;
//		n = nvA;
//		for (i = 0; i < nvB; i++) {
//			// ���B�еĽ��û����A�еĽ����ԣ�����������B�еĽ�����µ������еĽ����
//			if (newNodesB[i] == -1) {
//				newNodesB[i] = n;
//				n++;
//			} else {
//				// ����Ѿ���ԣ����������ԵĽ��ĸ���
//				nConnected++;
//			}
//
//		}
//
//		nBd = 0;
//		n = nb;
//		for (i = 0; i < meshB.nb; i++) {
//			if (newBdB[i] == -1) {
//				newBdB[i] = n;
//				n++;
//			} else {
//				nBd++;
//			}
//		}
//
//		nv = nvA + nvB - nConnected;
//		nt = nt + meshB.nt;
//		nb = nb + meshB.nb - nBd;
//
//		for (i = 0; i < nvB; i++) {
//			if (newNodesB[i] > nvA - 1) {
//				Node node = meshB.nodes.get(i);
//				nodes.add(new Node3D((Node3D) node));
//			}
//		}
//
//		for (i = 0; i < meshB.nt; i++) {
//			Element ele = meshB.elements.get(i);
//			int[] newEle = new int[3];
//			for (int kk = 0; kk < 3; kk++) {
//				newEle[kk] = newNodesB[ele.idx[kk]];
//			}
//			elements.add(new Tetrahedron4(newEle, ele.label));
//		}
//
//		for (i = 0; i < boundarys.size(); i++) {
//			// System.out.println(i + "\t" + nb + "\t" + boundarys.size());
//			if (ABd[i] == 1) {
//				boundarys.get(i).label = 0;
//			}
//		}
//
//		for (i = 0; i < meshB.nb; i++) {
//			if (newBdB[i] != 1) {
//				Element li = meshB.boundarys.get(i);
//				int[] newEle = new int[2];
//				for (int kk = 0; kk < 2; kk++) {
//					newEle[kk] = newNodesB[li.idx[kk]];
//				}
//				boundarys.add(new Line2(newEle, li.label));
//			}
//		}
//		return this;
//	}

//	public Mesh2DT3 add(Mesh2DT3... mesh) {
//		for (int i = 0; i < mesh.length; i++) {
//			this.add(mesh[i]);
//		}
//		return this;
//	}

	@Override
	public double[] getinteriorPointInElement(double[][] coord) {
		 double[] x = new double[3];
		 double[] y = new double[3];
		 for(int i = 0; i < 3; i++) {
			 x[i] = coord[i][0];
			 y[i] = coord[i][1];
		 }
		 return new double[]{ MVO.average(x), MVO.average(y) };

//		double[] aveCoord = new double[nDim];
//		for (int i = 0; i < nDim; i++) {
//			for (int j = 0; j < nPerElement; j++) {
//				aveCoord[i] += coord[j][i];
//			}
//			aveCoord[i] /= nPerElement;
//		}
//		return aveCoord;
	}

	@Override
	public double[] getinteriorPointOnBoundary(double[][] coord) {
		// double[] x = new double[2];
		// double[] y = new double[2];
		// for (int i = 0; i < 2; i++) {
		// x[i] = coord[i][0];
		// y[i] = coord[i][1];
		// }
		// return new double[]{ 0.5 * MVO.sum(x), 0.5 * MVO.sum(y)};

		double[] aveCoord = new double[nDim];
		for (int i = 0; i < nDim; i++) {
			for (int j = 0; j < nPerBoundary; j++) {
				aveCoord[i] += coord[j][i];
			}
			aveCoord[i] /= nPerBoundary;
		}
		return aveCoord;
	}

	@Override
	public int[][] getBoundarysFromElement(int i) {
		int[][] boundary = new int[nBoundaryPerElement][];
		int[] ele = elements[i];
		boundary[0] = new int[] { ele[0], ele[1] };
		boundary[1] = new int[] { ele[1], ele[2] };
		boundary[2] = new int[] { ele[2], ele[0] };
		return boundary;
	}

	/**
	 * �õ�����������Ԫ������ÿ��������Ӧ�ĵ�Ԫ��ţ�ʹ��������ϡ�����mat�洢��mat�е�i�е�j��Ԫ��val
	 * ��ʾ����Ԫ������������i��j����ɵı����ڵ�val����Ԫ��ע�⵽��ʱ�ߵı�����з���ģ�i-j�߶�Ӧ �ĵ�Ԫ��j-i�߶�Ӧ�ĵ�Ԫ�ǲ�ͬ�ġ�
	 * 
	 * @return ������ϡ�����
	 */
	public SIMatrix edgeFromElement() {
		SIMatrix edgeMat = new SIMatrix(nv);
		for (int i = 0; i < nt; i++) {
			int[] idx = elements[i];
			edgeMat.setElement(idx[0], idx[1], i);
			edgeMat.setElement(idx[1], idx[2], i);
			edgeMat.setElement(idx[2], idx[0], i);
		}
		return edgeMat;
	}

	public Mesh2T3 readMesh(String fileName, MeshFormat mf) {
		Mesh2T3 mesh = new Mesh2T3();
		Scanner sc = null;
		String line = "";
		int nv, nt, nb;
		try {
			sc = new Scanner(new File(fileName));
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}

		mesh.nDim = 2;
		mesh.nPerElement = 3;
		mesh.nPerBoundary = 2;
		mesh.nBoundaryPerElement = 3;

		switch (mf) {
			case msh:
				nv = sc.nextInt();
				nt = sc.nextInt();
				nb = sc.nextInt();
				mesh.initMesh0(nv, nt, nb);
				for (int i = 0; i < nv; i++) {
					double x = sc.nextDouble();
					double y = sc.nextDouble();
					mesh.nodes[i] = new double[] {x, y};
					mesh.nodeLabel[i] = sc.nextInt();
				}

				for (int i = 0; i < nt; i++) {
					int ti = sc.nextInt() - 1;
					int tj = sc.nextInt() - 1;
					int tk = sc.nextInt() - 1;
					mesh.elements[i] = new int[] {ti, tj, tk};
					mesh.elementLabel[i]  = sc.nextInt();
				}

				for (int i = 0; i < nb; i++) {
					int bi = sc.nextInt() - 1;
					int bj = sc.nextInt() - 1;
					mesh.boundaries[i] = new int[] {bi, bj};
					mesh.boundaryLabel[i] = sc.nextInt();
				}

				sc.close();
				break;
			case mmsh:
				nv = sc.nextInt();
				nt = sc.nextInt();
				mesh.nv = nv;
				mesh.nt = nt;
				mesh.nodes = new double[nv][];
				mesh.nodeLabel = new int[nv];
				mesh.elements = new int[nt][];
				mesh.elementLabel = new int[nt];
				for (int i = 0; i < nv; i++) {
					double x = sc.nextDouble();
					double y = sc.nextDouble();
					mesh.nodes[i] = new double[] {x, y};
					mesh.nodeLabel[i] = sc.nextInt();
				}

				for (int i = 0; i < nt; i++) {
					int ti = sc.nextInt() - 1;
					int tj = sc.nextInt() - 1;
					int tk = sc.nextInt() - 1;
					mesh.elements[i] = new int[] {ti, tj, tk};
					mesh.elementLabel[i]  = sc.nextInt();
				}
				mesh.extractBoundarysFromElements();
				sc.close();
				break;
			case Inp:
				do {
					line = sc.nextLine();
				}while(!line.startsWith("*Node"));
				line = sc.nextLine();
				ArrayList<double[] > nodes = new ArrayList<>();
				while (!line.startsWith("*Element")){
					String[] nodesCoord = line.split(",");
					double x = Double.parseDouble(nodesCoord[1].trim());
					double y = Double.parseDouble(nodesCoord[2].trim());
					nodes.add(new double[]{x, y});
					line = sc.nextLine();
				}
				nv = nodes.size();
				mesh.nv = nv;
				mesh.nodes = new double[nv][];
				mesh.nodeLabel = new int[nv];
				int idx = 0;
				for(double[] xy: nodes) {
					mesh.nodes[idx++] = xy.clone();
				}
				line = sc.nextLine();
				ArrayList<int[]> elements = new ArrayList<>();
				while(!line.startsWith("*Nset")) {
					String[] connectivity = line.split(",");
					int t1 = Integer.parseInt(connectivity[1].trim()) - 1;
					int t2 = Integer.parseInt(connectivity[2].trim()) - 1;
					int t3 = Integer.parseInt(connectivity[3].trim()) - 1;
					elements.add(new int[]{t1, t2, t3});
					line = sc.nextLine();
				}
				nt = elements.size();
				mesh.elements = new int[nt][];
				idx = 0;
				for(int[] ele: elements) {
					mesh.elements[idx++] = ele.clone();
				}
				mesh.elementLabel = new int[nt];
				mesh.extractBoundarysFromElements();
				//System.out.println(mesh.nv + "\t" + mesh.nt + "\t" + mesh.nb);
			default:
				break;
		}
		return mesh;
	}

	public Mesh2T3 readMesh(String xyFile, String eleFile) {
		Mesh2T3 mesh = new Mesh2T3();
		Scanner scXY = null;
		Scanner scEle = null;
		FileWriter fwOther = null;
		String line = "";
		try {
			scXY = new Scanner(new File(xyFile));
			scEle = new Scanner(new File(eleFile));
		} catch (IOException e) {
			e.printStackTrace();
		}

		ArrayList<double[]> nodes = new ArrayList<>();
		while(scXY.hasNext()) {
			line = scXY.nextLine();
			String[] xy = line.split(",");
			double t1 = Double.parseDouble(xy[0]);
			double t2 = Double.parseDouble(xy[1]);
			nodes.add(new double[] {t1, t2});
		}
		ArrayList<int[]> eles = new ArrayList<>();
		while(scEle.hasNext()) {
			line = scEle.nextLine();
			String[] ele = line.split(",");
			int t1 = Integer.parseInt(ele[0]);
			int t2 = Integer.parseInt(ele[1]);
			int t3 = Integer.parseInt(ele[2]);
			eles.add(new int[] {t1, t2, t3});
		}
		int nv = nodes.size();
		int nt = eles.size();
		int nb = 0;
		mesh.initMesh0(nv, nt, nb);
		for(int i = 0; i < nv; i++) {
			mesh.nodes[i] = nodes.get(i);
		}
		for(int i = 0; i < nt; i++) {
			mesh.elements[i] = eles.get(i);
		}
		mesh.extractBoundarysFromElements();

		return mesh;
	}

	public Mesh2T3 square2D(int xNum, int yNum, T3Type type) {
		int nv = (xNum + 1) * (yNum + 1);
		int nt = 2 * xNum * yNum;
		int nb = 2 * (xNum + yNum);
		nv = nv;
		nt = nt;
		nb = nb;

		initMesh0(nv, nt, nb);
		int count = 0;
		for (int i = 0; i < yNum + 1; i++) {
			for (int j = 0; j < xNum + 1; j++) {
				nodes[count++] = new double[] {(j + 0.0) / xNum,(i + 0.0) / yNum};
			}
		}

		createSquareElements(xNum, yNum, type);
		createSquareBoundarys(xNum, yNum, Direct.X);
		return this;
	}

	public Mesh2T3 square2D(int xNum, int yNum, T3Type type, Direct direct) {
		if (direct == Direct.X) {
			return square2D(xNum, yNum, type);
		} else if (direct == Direct.Y) {
			Mesh2T3 mesh = new Mesh2T3();
			int nv = (xNum + 1) * (yNum + 1);
			int nt = 2 * xNum * yNum;
			int nb = 2 * (xNum + yNum);
			mesh.nv = nv;
			mesh.nt = nt;
			mesh.nb = nb;
			int count = 0;
			mesh.initMesh0(nv, nt, nb);
			for (int j = 0; j < xNum + 1; j++) {
				for (int i = 0; i < yNum + 1; i++) {
					mesh.nodes[count++] = new double[] {(j + 0.0) / xNum, (i + 0.0) / yNum};
				}
			}
			createSquareElements(yNum, xNum, type, direct);
			createSquareBoundarys(yNum, xNum, direct);
			return mesh;
		} else {
			throw new IllegalArgumentException("The parameter Direct must only be Direct.X or Direct.Y!");
		}
	}

	public Mesh2T3 square(int xNum, int yNum) {
		int t1, t2, t3, t4, t5;
		nv = 2 * xNum * yNum + xNum + yNum + 1;
		nt = 4 * xNum * yNum;
		nb = 2 * (xNum + yNum);
		initMesh0(nv, nt, nb);
		int count = 0;
		for(int j = 0; j < yNum; j++) {
			for(int i = 0; i < xNum + 1; i++) {
				nodes[count++] = new double[] {(i + 0.0) / xNum, (j + 0.0) / yNum};
			}
			for(int i = 0; i < xNum; i++) {
				nodes[count++] = new double[] {(i + 0.5) / xNum, (j + 0.5) / yNum};
			}
		}

		for(int i = 0; i < xNum + 1; i++) {
			nodes[count++] = new double[]{(i + 0.0) / xNum, 1.0};
		}

		count = 0;
		for(int j = 0; j < yNum; j++) {
			for(int i = 0; i < xNum; i++) {
				t1 = i + j * (2 * xNum + 1);
				t2 = t1 + 1;
				t3 = t2 + 2 * xNum + 1;
				t4 = t3 - 1;
				t5 = t1 + xNum + 1;
				elements[count++] = new int[] {t1, t2, t5};
				elements[count++] = new int[] {t2, t3, t5};
				elements[count++] = new int[] {t3, t4, t5};
				elements[count++] = new int[] {t4, t1, t5};
			}
		}

		count  = 0;
		for(int i = 0; i < xNum; i++) {
			boundaries[count] = new int[] {i, i + 1};
			boundaryLabel[count++] = 1;
		}

		int local = xNum;
		for(int i = 0; i < yNum; i++) {
			boundaries[count] = new int[] {local, local + 2 * xNum + 1};
			boundaryLabel[count++] = 2;
			local += 2 * xNum + 1;
		}

		local = 2 * xNum * yNum + xNum + yNum;
		for(int i = 0; i < xNum; i++) {
			boundaries[count] = new int[] {local, local - 1};
			boundaryLabel[count++] = 3;
			local--;
		}

		local = yNum * (2 * xNum + 1);
		for(int i = 0; i < yNum; i++) {
			boundaries[count] = new int[] {local, local - 2 * xNum - 1};
			boundaryLabel[count++] = 4;
			local -= 2 * xNum + 1;
		}
		return this;
	}

	public Mesh2T3 quadrilateral(double[] QuadX, double[] QuadY, int xNum, int yNum) {
		int n, i, j, xNodes, yNodes, ti, tj, tk, tl;
		double xi, et, x, y;
		double[] h = new double[4];
		double[] dXi = new double[4];
		double[] dEt = new double[4];
		double[] hh = new double[8];
		double[] ddXi = new double[8];
		double[] ddEt = new double[8];
		double[] tji = new double[2];
		double[] tli = new double[2];
		n = QuadX.length;
		xNodes = xNum + 1;
		yNodes = yNum + 1;
		int nv = xNodes * yNodes;
		int nt = 2 * xNum * yNum;
		int nb = 2 * (xNum + yNum);
		initMesh0(nv, nt, nb);
		int count = 0;
		if (n == 4) {
			for (j = 0; j < yNodes; j++) {
				for (i = 0; i < xNodes; i++) {
					xi = 2 * (i + 0.0) / xNum - 1.0;
					et = 2 * (j + 0.0) / yNum - 1.0;
					CG.Q4ShapeFunction(xi, et, h, dXi, dEt);
					x = MVO.dot_product(h, QuadX);
					y = MVO.dot_product(h, QuadY);
					nodes[count++] = new double[] {x, y};
				}
			}
		} else if (n == 8) {
			for (j = 0; j < yNodes; j++) {
				for (i = 0; i < xNodes; i++) {
					xi = 2 * (i + 0.0) / xNum - 1.0;
					et = 2 * (j + 0.0) / yNum - 1.0;
					CG.Q8ShapeFunction(xi, et, hh, ddXi, ddEt);
					x = MVO.dot_product(hh, QuadX);
					y = MVO.dot_product(hh, QuadY);
					nodes[count++] = new double[] {x, y};
				}
			}
		} else
			throw new IllegalArgumentException("The length of the array must be 4 or 8!");

		count = 0;
		for (j = 0; j < yNum; j++) {
			for (i = 0; i < xNum; i++) {
				ti = i + j * xNodes;
				tj = ti + 1;
				tk = tj + xNodes;
				tl = tk - 1;
				for (int kk = 0; kk < 2; kk++) {
					tji[kk] = nodes[tj][kk] - nodes[ti][kk];
					tli[kk] = nodes[tl][kk] - nodes[ti][kk];
				}
				if (MVO.dot_product(tji, tli) > 0.0) {
					elements[count++] = new int[] {ti, tj, tl};
					elements[count++] = new int[] {tj, tk, tl};
				} else {
					elements[count++] = new int[] {tj, tk, ti};
					elements[count++] = new int[] {tk, tl, ti};
				}
			}
		}

		createSquareBoundarys(xNum, yNum);
		return this;
	}

	public Mesh2T3 ring(double ra, double rb, int nR, int nL) {
		double[] theta = MVO.linspace(0, 2 * PI, nR + 1);
		nv = (nL + 1) * nR;
		nt = 2 * nR * nL;
		nb = 2 * nR;
		initMesh0(nv, nt, nb);
		double dr = (rb - ra) / nL;
		double x, y;
		int count = 0;
		for (int i = 0; i < nL + 1; i++) {
			for (int j = 0; j < nR; j++) {
				x = (ra + i * dr) * cos(theta[j]);
				y = (ra + i * dr) * sin(theta[j]);
				// System.out.println(MVector.toString(x, y));
				nodes[count++] = new double[] {x, y};
			}
		}

		count = 0;
		int L1, L2, L3, L4;
		for (int i = 0; i < nL; i++) {
			for (int j = 0; j < nR - 1; j++) {
				L1 = j + i * nR;
				L2 = L1 + nR;
				L3 = L2 + 1;
				L4 = L1 + 1;
				elements[count++] = new int[] {L1, L2, L3};
				elements[count++] = new int[] {L1, L3, L4};
			}
			L1 = i * nR + nR - 1;
			L2 = L1 + nR;
			L3 = (i + 1) * nR;
			L4 = i * nR;
			elements[count++] = new int[] {L1, L2, L3};
			elements[count++] = new int[] {L1, L3, L4};
		}

		count = 0;
		for (int i = 0; i < nR - 1; i++) {
			boundaries[count] = new int[] {i + 1, i};
			boundaryLabel[count++] = 1;
		}
		boundaries[count] = new int[] {0, nR - 1};
		boundaryLabel[count++] = 1;

		int local = nL * nR;
		for (int i = 0; i < nR - 1; i++) {
			boundaries[count] = new int[] {local, local + 1};
			boundaryLabel[count++] = 2;
			local++;
		}

		boundaries[count] = new int[] {nR * (nL + 1) - 1, nL * nR};
		boundaryLabel[count] = 2;
		return this;
	}

	public Mesh2T3 triangle(double[] x, double[] y, int num) {
		double avX = MVO.average(x);
		double avY = MVO.average(y);
		double[] quadX1 = {x[0], 0.5 * (x[0] + x[1]), avX, 0.5 * (x[0] + x[2])};
		double[] quadY1 = {y[0], 0.5 * (y[0] + y[1]), avY, 0.5 * (y[0] + y[2])};
		quadrilateral(quadX1, quadY1, num / 2, num / 2);

		double[] quadX2 = {x[1], 0.5 * (x[1] + x[2]), avX, 0.5 * (x[0] + x[1])};
		double[] quadY2 = {y[1], 0.5 * (y[1] + y[2]), avY, 0.5 * (y[0] + y[1])};
		Mesh2T3 mesh2 = new Mesh2T3();
		mesh2.quadrilateral(quadX2, quadY2, num/2, num/2);
		this.add(mesh2);

		double[] quadX3 = {x[2], 0.5 * (x[2] + x[0]), avX, 0.5 * (x[1] + x[2])};
		double[] quadY3 = {y[2], 0.5 * (y[2] + y[0]), avY, 0.5 * (y[1] + y[2])};
		Mesh2T3 mesh3 = new Mesh2T3();
		mesh3.quadrilateral(quadX3, quadY3, num/2, num/2);
		this.add(mesh3);
		return this;
	}

	public Mesh2T3 fan(double ra, double rb, double rc, double rd, double theta1, double theta2,
					   int rNum, int thetaNum, T3Type type) {
		int nv = (rNum + 1) * (thetaNum + 1);
		int nt = 2 * rNum * thetaNum;
		int nb = 2 * (rNum + thetaNum);
		initMesh0(nv, nt, nb);
		double rab = rb - ra;
		double rdc = rc - rd;
		double dr0 = rab / rNum;
		double dr1 = rdc / rNum;

		int count = 0;
		for (int i = 0; i < thetaNum + 1; i++) {
			double dr = dr0 + i * (dr1 - dr0) / thetaNum;
			double R0 = ra + i * (rd - ra) / thetaNum;
			double theta = theta1 + i * (theta2 - theta1) / thetaNum;
			for (int j = 0; j < rNum + 1; j++) {
				double R = R0 + dr * j;
				nodes[count++] = new double[] { R * cos(theta), R * sin(theta)};
			}
		}

		createSquareElements(rNum, thetaNum, type);
		createSquareBoundarys(rNum, thetaNum, Direct.X);
		return this;
	}

	/**
	 * 锟皆撅拷锟叫撅拷锟斤拷锟斤拷锟剿结构锟斤拷锟斤拷锟斤拷悖拷锟斤拷锟斤拷涞ピ拷锟斤拷锟斤拷锟斤拷锟斤拷呓锟斤拷锟斤拷锟斤拷锟�
	 *
	 *            要锟斤拷锟缴碉拷锟斤拷锟斤拷
	 * @param xNum
	 *            x锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟�
	 * @param yNum
	 *            y锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟�
	 * @param type
	 *            锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟轿凤拷锟斤拷锟斤拷锟斤拷
	 */
	private void createSquareElements(int xNum, int yNum, T3Type type) {
		int t1, t2, t3, t4;
		int count = 0;
		for (int i = 0; i < yNum; i++) {
			for (int j = 0; j < xNum; j++) {
				t1 = i * (xNum + 1) + j;
				t2 = t1 + 1;
				t3 = t1 + xNum + 1;
				t4 = t3 + 1;
				switch (type) {
					case Left:
						elements[count++] = new int[] {t1, t2, t4};
						elements[count++] = new int[] {t1, t4, t3};
						break;
					case Right:
						elements[count++] = new int[] {t1, t2, t3};
						elements[count++] = new int[] {t2, t4, t3};
						break;
					case Cross1:
						if ((i + j) % 2 == 0) {
							elements[count++] = new int[] {t1, t2, t4};
							elements[count++] = new int[] {t1, t4, t3};
						} else {
							elements[count++] = new int[] {t1, t2, t3};
							elements[count++] = new int[] {t2, t4, t3};
						}
						break;
					case Cross2:
						if ((i + j) % 2 == 1) {
							elements[count++] = new int[] {t1, t2, t4};
							elements[count++] = new int[] {t1, t4, t3};
						} else {
							elements[count++] = new int[] {t1, t2, t3};
							elements[count++] = new int[] {t2, t4, t3};
						}
						break;
					default:
						throw new IllegalArgumentException("Unexpected T3Type!");
				}
			}
		}
	}

	private void createSquareElements(int xNum, int yNum, T3Type type, Direct direct) {
		if (direct == Direct.X) {
			createSquareElements(xNum, yNum, type);
		} else if (direct == Direct.Y) {
			int count = 0;
			for (int i = 0; i < yNum; i++) {
				for (int j = 0; j < xNum; j++) {
					int ti1, tj1, tk1, ti2, tj2, tk2;
					switch (type) {
						case Left:
							ti1 = i * (xNum + 1) + j;
							tj1 = ti1 + 1;
							tk1 = tj1 + xNum + 1;
							ti2 = i * (xNum + 1) + j;
							tj2 = ti2 + xNum + 1 + 1;
							tk2 = tj2 - 1;
							break;
						case Right:
							ti1 = i * (xNum + 1) + j;
							tj1 = ti1 + 1;
							tk1 = ti1 + xNum + 1;
							ti2 = i * (xNum + 1) + j + 1;
							tj2 = ti2 + xNum + 1;
							tk2 = tj2 - 1;
							break;
						case Cross1:
							if ((i + j) % 2 == 0) {
								ti1 = i * (xNum + 1) + j;
								tj1 = ti1 + 1;
								tk1 = ti1 + xNum + 1;
								ti2 = i * (xNum + 1) + j + 1;
								tj2 = ti2 + xNum + 1;
								tk2 = tj2 - 1;
							} else {
								ti1 = i * (xNum + 1) + j;
								tj1 = ti1 + 1;
								tk1 = tj1 + xNum + 1;
								ti2 = i * (xNum + 1) + j;
								tj2 = ti2 + xNum + 1 + 1;
								tk2 = tj2 - 1;
							}
							break;
						case Cross2:
							if ((i + j) % 2 == 1) {
								ti1 = i * (xNum + 1) + j;
								tj1 = ti1 + 1;
								tk1 = ti1 + xNum + 1;
								ti2 = i * (xNum + 1) + j + 1;
								tj2 = ti2 + xNum + 1;
								tk2 = tj2 - 1;
							} else {
								ti1 = i * (xNum + 1) + j;
								tj1 = ti1 + 1;
								tk1 = tj1 + xNum + 1;
								ti2 = i * (xNum + 1) + j;
								tj2 = ti2 + xNum + 1 + 1;
								tk2 = tj2 - 1;
							}
							break;
						default:
							throw new IllegalArgumentException("Unexpected T3Type!");
					}
					elements[count++] = new int[] {ti1, tk1, tj1};
					elements[count++] = new int[] {ti2, tk2, tj2};
				}
			}
		} else {
			throw new IllegalArgumentException("The parameter Direct must only be Direct.X or Direct.Y!");
		}
	}

	private void createSquareBoundarys(int xNum, int yNum) {
		int count = 0;
		for (int i = 0; i < xNum; i++) {
			boundaries[count] = new int[] {i, i + 1};
			boundaryLabel[count++] = 1;
		}
		int local = xNum;
		for (int i = 0; i < yNum; i++) {
			boundaries[count] = new int[] {local, local + xNum + 1};
			boundaryLabel[count++] = 2;
			local += xNum + 1;
		}

		local = xNum + yNum * (xNum + 1);
		for (int i = 0; i < xNum; i++) {
			boundaries[count] = new int[] {local, local - 1};
			boundaryLabel[count++] = 3;
			local--;
		}
		local = yNum * (xNum + 1);
		for (int i = 0; i < yNum; i++) {
			boundaries[count] = new int[] {local, local - (xNum + 1)};
			boundaryLabel[count++] = 4;
			local -= xNum + 1;
		}
	}

	private void createSquareBoundarys(int xNum, int yNum, Direct direct) {
		if (direct == Direct.X) {
			createSquareBoundarys(xNum, yNum);
		} else if (direct == Direct.Y) {
			int count = 0;
			for (int i = 0; i < xNum; i++) {
				boundaries[count] = new int[] {i + 1, i};
				boundaryLabel[count++] = 4;
			}
			int local = xNum;
			for (int i = 0; i < yNum; i++) {
				boundaries[count] = new int[] {local + xNum + 1, local};
				boundaryLabel[count++] = 3;
				local += xNum + 1;
			}

			local = xNum + yNum * (xNum + 1);
			for (int i = 0; i < xNum; i++) {
				boundaries[count] = new int[] {local - 1, local};
				boundaryLabel[count++] = 2;
				local--;
			}

			local = yNum * (xNum + 1);
			for (int i = 0; i < yNum; i++) {
				boundaries[count] = new int[] {local - (xNum + 1), local};
				boundaryLabel[count++] = 1;
				local -= xNum + 1;
			}
		} else {
			throw new IllegalArgumentException("The parameter Direct must only be Direct.X or Direct.Y!");
		}

	}

	public Mesh2T3 fromRectangleNodes(double[][] nodes, int xNum, int yNum, T3Type type) {
		int nv = nodes.length;
		ParamCheck.checkEqual(nv, (xNum + 1) * (yNum + 1));
		int nt = 2 * xNum * yNum;
		int nb = 2 * (xNum + yNum);
		initMesh0(nv, nt, nb);
		for(int i = 0; i < nv; i++) {
			nodes[i] = nodes[i].clone();
		}

		createSquareElements(xNum, yNum, type);
		createSquareBoundarys(xNum, yNum, xNum <= yNum? Direct.X : Direct.Y);
		return this;
	}

	public static Mesh2T3 triangulate(double[][] pxy) {
		int nv = pxy.length - 3;
		int[][] edges = null;
		boolean complete[] = null;//用来存储三角形还是否需要被考虑外接圆
		int nedge = 0;
		int trimax, emax = 200;
		int status;
		int[][]  triangles = null;
		boolean inside;
		double 	xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
		double 	xmin, xmax, ymin, ymax, xmid, ymid;
		double 	dx, dy, dmax;
		int ntri = 0;

		//给complete列表分配内存
		trimax = 4 * nv;
		complete = new boolean[trimax];
		for (int i = 0; i < trimax; i++)  {
			complete[i] = false;
		}

		//给边缘列表分配内存
		emax = 4 * nv;
		edges = new int[emax][];
		for (int i =0; i < emax; i++) {
			edges[i] = new int[2];
		}

		/*
		    找到所有点横纵坐标的最大最小值
		*/
		xmax = xmin = pxy[0][0];
		ymax = ymin = pxy[0][1];
		for (int i = 1; i < nv; i++)
		{
			if (pxy[i][0] < xmin) {
				xmin = pxy[i][0];
			}
			if (pxy[i][0] > xmax) {
				xmax = pxy[i][0];
			}
			if (pxy[i][1] < ymin) {
				ymin = pxy[i][1];
			}
			if (pxy[i][1] > ymax) {
				ymax = pxy[i][1];
			}
		}
		dx = xmax - xmin;
		dy = ymax - ymin;
		dmax = (dx > dy) ? dx : dy;
		xmid = 0.5 * (xmax + xmin);
		ymid = 0.5 * (ymax + ymin);

		/* 设定超三角形（一个包含所有样本点的大三角形）
		   将超三角形的顶点加入到顶点列表pxy的最后
		   超三角形是三角形列表中的第一个三角形
		*/

		pxy[nv][0] = xmid - 2.0 * dmax;
		pxy[nv][1] = ymid - dmax;
		pxy[nv + 1][0] = xmid;
		pxy[nv + 1][1] = ymid + 2.0 * dmax;
		pxy[nv + 2][0] = xmid + 2.0 * dmax;
		pxy[nv + 2][1] = ymid - dmax;

		int[][]	v = new int [nv*3][];
		v[0] = new int[]{nv, nv + 1, nv + 2};
		complete[0] = false;
		ntri = 1;

		/*
			每次考虑一个点加入到已有的网格中去
		*/
		for (int i = 0; i < nv; i++) {
			xp = pxy[i][0];
			yp = pxy[i][1];
			nedge = 0;

			/*
				Set up the edge buffer.
				如果点(xp,yp)在这个三角形的外接圆内，则把这个三角形的三条边加入到edge buffer中去
                然后把这个三角形从三角形列表中删去
			*/
			double[] circle = new double[3];
			for (int j = 0;j < ntri; j++) {
				if (complete[j])
					continue;
				int[] tri = v[j];
				double[] x = {pxy[tri[0]][0], pxy[tri[1]][0], pxy[tri[2]][0]};
				double[] y = {pxy[tri[0]][1], pxy[tri[1]][1], pxy[tri[2]][1]};
				double[] xyp = {xp, yp};
				inside = CG.CircumCircle(xyp, x, y, circle);
				xc = circle[0];
				yc = circle[1];
				r = circle[2];

				//如果圆上横坐标最大的点的横坐标都比当前xp小，那么这个三角形以后将不用考虑
				if (xc + r < xp) {
					complete[j] = true;
				}

				if (inside) {
					//检查是否超出边缘列表，如果超出了就给列表长度加100
					if (nedge + 3 >= emax) {
						emax += 100;
						int[][] edges_n = new int[emax][];
						for (int ie = 0; ie < emax; ie++) {
							edges_n[ie] = new int[2];
						}
						System.arraycopy(edges, 0, edges_n, 0, edges.length);
						edges = edges_n;
					}
					//把该三角形存入边缘列表
					edges[nedge + 0][0] = v[j][0];
					edges[nedge + 0][1] = v[j][1];

					edges[nedge + 1][0] = v[j][1];
					edges[nedge + 1][1] = v[j][2];

					edges[nedge + 2][0] = v[j][2];
					edges[nedge + 2][1] = v[j][0];
					nedge += 3;  //边界条数加3
					v[j][0] = v[ntri - 1][0];
					v[j][1] = v[ntri - 1][1];
					v[j][2] = v[ntri - 1][2];
					complete[j] = complete[ntri - 1];
					ntri--;
					j--;
				}
			}

			/*
				Tag multiple edges
				Note: if all triangles are specified anticlockwise then all
				interior edges are opposite pointing in direction.
			*/
			for (int j = 0;j < nedge - 1;j++)
			{
				for (int k = j + 1;k < nedge; k++)
				{
					//对比边是否一样，一样的删去（赋值-1）
					if ((edges[j][0] == edges[k][1]) && (edges[j][1] == edges[k][0])) {
						edges[j][0] = -1;
						edges[j][1] = -1;
						edges[k][0] = -1;
						edges[k][1] = -1;
					}
					if ((edges[j][0] == edges[k][0]) && (edges[j][1] == edges[k][1])) {
						edges[j][0] = -1;
						edges[j][1] = -1;
						edges[k][0] = -1;
						edges[k][1] = -1;
					}
				}
			}

			int nb = 0;
			System.out.println("Edges");
			/*
				在当前点形成新的三角形
				跳过重合的边
				所有边都是顺时针的
			*/
			for (int j = 0; j < nedge; j++)
			{
				System.out.println(Arrays.toString(edges[j]));
				if (edges[j][0] == -1 || edges[j][1] == -1)
					continue;
				if (ntri >= trimax) {
					return null;
				}
				v[ntri][0] = edges[j][0];
				v[ntri][1] = edges[j][1];
				v[ntri][2] = i;
				complete[ntri] = false;
				ntri ++;
			}
		}


		/*
			删去和超三角形有公共顶点的三角形
			是那些顶点编号大于等于nv的三角形
		*/
		for (int i = 0; i < ntri;i ++) {
			if (v[i][0] >= nv || v[i][1] >= nv || v[i][2] >= nv) {
				v[i] = v[ntri - 1];
				ntri--;
				i--;
			}
		}
		Mesh2T3 mesh = new Mesh2T3();
		mesh.nv = nv;
		mesh.nt = ntri;
		mesh.nodes = new double[nv][];
		for(int i = 0; i < nv; i++) {
			mesh.nodes[i] = pxy[i].clone();
		}
		for(int i = 0; i < ntri; i++) {
			mesh.elements[i] = v[i].clone();
		}
		return mesh;
	}

//	public void adaptRefine(double[] errorIndicator, double ratio, ElementFunc bdFunc, double[] param) {
//		int oldnv = nv;
//		int oldNt = nt;
//		int oldnb = nb;
//		int[] refineIdx = new int[nt];
//		int[] eleIdx = new int[3];
//		ArrayList<Node> oldNodes = new ArrayList<>();
//		ArrayList<Element> oldElements = new ArrayList<>();
//		ArrayList<Element> oldBoundarys = new ArrayList<>();
//		nodes.stream().forEach(ele -> oldNodes.add(new Node2D((Node2D) ele)));
//		elements.stream().forEach(ele -> oldElements.add(new Triangle3((Triangle3) ele)));
//		boundarys.stream().forEach(ele -> oldBoundarys.add(new Line2((Line2) (ele))));
//
//		SIMatrix B = new SIMatrix(nv);
//		double maxError = MVO.max(errorIndicator);
//		for (int i = 0; i < nt; i++) {
//			int[] idx = getElement(i).idx;
//			if (errorIndicator[i] > maxError * ratio) {
//				B.setElement(idx[0], idx[1], 1);
//				B.setElement(idx[1], idx[0], 1);
//				B.setElement(idx[1], idx[2], 1);
//				B.setElement(idx[2], idx[1], 1);
//				B.setElement(idx[2], idx[0], 1);
//				B.setElement(idx[0], idx[2], 1);
//				refineIdx[i] = 3;
//			}
//		}
//		B.sort();
//
//		int pNNZ = 0;
//		double[][] coord = new double[3][2];
//		double[] length = new double[3];
//		Pair<Integer, Double> maxLocVal = new Pair<>();
//		Pair<Integer, Integer> minLocValI = new Pair<>();
//		Pair<Integer, Integer> maxLocValI = new Pair<>();
//		while (true) {
//			pNNZ = B.getNNZ();
//			for (int i = 0; i < nt; i++) {
//				Element ele = getElement(i);
//				int[] idx = ele.idx;
//				getCoordPerElement(ele, coord);
//				eleIdx[0] = B.getElement(idx[0], idx[1]);
//				eleIdx[1] = B.getElement(idx[1], idx[2]);
//				eleIdx[2] = B.getElement(idx[2], idx[0]);
//				CG.length2DT3(coord, length);
//				MVO.max(length, maxLocVal);
//				int maxIdx = maxLocVal.getFirst();
//				MVO.min(eleIdx, minLocValI);
//				MVO.max(eleIdx, maxLocValI);
//				int sumEleIdx = MVO.sum(eleIdx[0]);
//				if (sumEleIdx == 2) {
//					int zeroIdx = minLocValI.getFirst();
//					if (maxIdx == zeroIdx) {
//						refineIdx[i] = 3;
//						if (maxIdx == 0) {
//							B.setElement(idx[0], idx[1], 1);
//							B.setElement(idx[1], idx[0], 1);
//						} else if (maxIdx == 1) {
//							B.setElement(idx[1], idx[2], 1);
//							B.setElement(idx[2], idx[1], 1);
//						} else {
//							B.setElement(idx[2], idx[0], 1);
//							B.setElement(idx[0], idx[2], 1);
//						}
//					}
//				} else if (sumEleIdx == 1) {
//					int zeroIdx = maxLocValI.getFirst();
//					if (zeroIdx != maxIdx) {
//						refineIdx[i] = 2;
//						if (maxIdx == 0) {
//							B.setElement(idx[0], idx[1], 1);
//							B.setElement(idx[1], idx[0], 1);
//						} else if (maxIdx == 1) {
//							B.setElement(idx[1], idx[2], 1);
//							B.setElement(idx[2], idx[1], 1);
//						} else {
//							B.setElement(idx[2], idx[0], 1);
//							B.setElement(idx[0], idx[2], 1);
//						}
//					}
//				}
//			}
//			if (pNNZ == B.getNNZ()) {
//				break;
//			}
//		}
//
//		for (int i = 0; i < nt; i++) {
//			int[] idx = elements.get(i).idx;
//			eleIdx[0] = B.getElement(idx[0], idx[1]);
//			eleIdx[1] = B.getElement(idx[1], idx[2]);
//			eleIdx[2] = B.getElement(idx[2], idx[0]);
//			refineIdx[i] = MVO.sum(eleIdx);
//		}
//
//		int newNv = nv + B.getNNZ() / 2;
//		int newNt = nt + MVO.sum(refineIdx);
//		nv = newNv;
//		nt = newNt;
//		int nIdx = nv - 1;
//		int eIdx = 0;
//		int[] insertEleIdx = new int[3];
//		double insertX, insertY;
//		elements.clear();
//		for (int i = 0; i < oldNt; i++) {
//			Element oldEle = oldElements.get(i);
//			// getCoordPerElement(oldEle, coord);
//			for (int j = 0; j < 3; j++) {
//				coord[j] = oldNodes.get(oldEle.idx[j]).coord.clone();
//			}
//			int[] idx = oldEle.idx;
//			eleIdx[0] = B.getElement(idx[0], idx[1]);
//			eleIdx[1] = B.getElement(idx[1], idx[2]);
//			eleIdx[2] = B.getElement(idx[2], idx[0]);
//			if (refineIdx[i] == 3) {
//				insertEleIdx = eleIdx.clone();
//				if (insertEleIdx[0] == 1) {
//					nIdx++;
//					insertX = 0.5 * (coord[0][0] + coord[1][0]);
//					insertY = 0.5 * (coord[0][1] + coord[1][1]);
//					nodes.add(new Node2D(insertX, insertY));
//					insertEleIdx[0] = nIdx;
//					B.setElement(idx[0], idx[1], nIdx);
//					B.setElement(idx[1], idx[0], nIdx);
//				}
//				if (insertEleIdx[1] == 1) {
//					nIdx++;
//					insertX = 0.5 * (coord[1][0] + coord[2][0]);
//					insertY = 0.5 * (coord[1][1] + coord[2][1]);
//					nodes.add(new Node2D(insertX, insertY));
//					insertEleIdx[1] = nIdx;
//					B.setElement(idx[1], idx[2], nIdx);
//					B.setElement(idx[2], idx[1], nIdx);
//				}
//
//				if (insertEleIdx[2] == 1) {
//					nIdx++;
//					insertX = 0.5 * (coord[2][0] + coord[0][0]);
//					insertY = 0.5 * (coord[2][1] + coord[0][1]);
//					nodes.add(new Node2D(insertX, insertY));
//					insertEleIdx[2] = nIdx;
//					B.setElement(idx[2], idx[0], nIdx);
//					B.setElement(idx[0], idx[2], nIdx);
//				}
//				elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[2], oldEle.label));
//				elements.add(new Triangle3(idx[1], insertEleIdx[1], insertEleIdx[0], oldEle.label));
//				elements.add(new Triangle3(idx[2], insertEleIdx[2], insertEleIdx[1], oldEle.label));
//				elements.add(new Triangle3(insertEleIdx, oldEle.label));
//			} else if (refineIdx[i] == 2) { // ���ﴦ�������������߲����ʱ�ļ���
//				CG.length2DT3(coord, length);
//				MVO.max(length, maxLocVal);
//				int maxIdx = maxLocVal.getFirst();
//				if (eleIdx[0] == 0) {
//					insertEleIdx[0] = eleIdx[1];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[1][0] + coord[2][0]);
//						insertY = 0.5 * (coord[1][1] + coord[2][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[0] = nIdx;
//						B.setElement(idx[2], idx[1], nIdx);
//						B.setElement(idx[1], idx[2], nIdx);
//					}
//					insertEleIdx[1] = eleIdx[2];
//					if (insertEleIdx[1] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[2][0] + coord[0][0]);
//						insertY = 0.5 * (coord[2][1] + coord[0][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[1] = nIdx;
//					}
//
//					if (maxIdx == 1) {
//						elements.add(new Triangle3(idx[0], idx[1], insertEleIdx[0], oldEle.label));
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[2], insertEleIdx[1], oldEle.label));
//					} else if (maxIdx == 2) {
//						elements.add(new Triangle3(idx[0], idx[1], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(idx[1], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[2], insertEleIdx[1], oldEle.label));
//					} else {
//						elements.add(new Triangle3(idx[0], idx[1], insertEleIdx[0], oldEle.label));
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[2], insertEleIdx[1], oldEle.label));
//					}
//				} else if (eleIdx[1] == 0) {
//					insertEleIdx[0] = eleIdx[0];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[0][0] + coord[1][0]);
//						insertY = 0.5 * (coord[0][1] + coord[1][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[0] = nIdx;
//						B.setElement(idx[0], idx[1], nIdx);
//						B.setElement(idx[1], idx[0], nIdx);
//					}
//					insertEleIdx[1] = eleIdx[2];
//					if (insertEleIdx[1] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[0][0] + coord[2][0]);
//						insertY = 0.5 * (coord[0][1] + coord[2][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[1] = nIdx;
//						B.setElement(idx[0], idx[2], nIdx);
//						B.setElement(idx[2], idx[0], nIdx);
//					}
//					if (maxIdx == 0) {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], idx[2], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[2], insertEleIdx[1], oldEle.label));
//					} else if (maxIdx == 2) {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], insertEleIdx[2], oldEle.label));
//						elements.add(new Triangle3(idx[1], idx[2], insertEleIdx[1], oldEle.label));
//					} else {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], idx[2], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[2], insertEleIdx[1], oldEle.label));
//					}
//				} else {
//					insertEleIdx[0] = eleIdx[0];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[0][0] + coord[1][0]);
//						insertY = 0.5 * (coord[0][1] + coord[1][1]);
//						insertEleIdx[0] = nIdx;
//						nodes.add(new Node2D(insertX, insertY));
//						B.setElement(idx[0], idx[1], nIdx);
//						B.setElement(idx[1], idx[0], nIdx);
//					}
//					insertEleIdx[1] = eleIdx[1];
//					if (insertEleIdx[1] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[1][0] + coord[2][0]);
//						insertY = 0.5 * (coord[1][1] + coord[2][1]);
//						insertEleIdx[1] = nIdx;
//						nodes.add(new Node2D(insertX, insertY));
//						B.setElement(idx[1], idx[2], nIdx);
//						B.setElement(idx[2], idx[1], nIdx);
//					}
//					if (maxIdx == 0) {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], idx[2], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], insertEleIdx[1], idx[2]));
//					} else if (maxIdx == 1) {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(idx[0], insertEleIdx[1], idx[2], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], insertEleIdx[1], oldEle.label));
//					} else {
//						elements.add(new Triangle3(idx[0], insertEleIdx[0], idx[2], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], idx[1], insertEleIdx[1], oldEle.label));
//						elements.add(new Triangle3(insertEleIdx[0], insertEleIdx[1], idx[2]));
//					}
//				}
//			} else if (refineIdx[i] == 1) { // ���ﴦ�������������߲����ʱ�ļ���
//				if (eleIdx[0] != 0) {
//					insertEleIdx[0] = eleIdx[0];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[0][0] + coord[1][0]);
//						insertY = 0.5 * (coord[0][1] + coord[1][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[0] = nIdx;
//						B.setElement(idx[0], idx[1], nIdx);
//						B.setElement(idx[1], idx[0], nIdx);
//					}
//					elements.add(new Triangle3(idx[0], insertEleIdx[0], idx[2], oldEle.label));
//					elements.add(new Triangle3(insertEleIdx[0], idx[1], idx[2], oldEle.label));
//				} else if (eleIdx[1] != 0) {
//					insertEleIdx[0] = eleIdx[1];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[1][0] + coord[2][0]);
//						insertY = 0.5 * (coord[1][1] + coord[2][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[0] = nIdx;
//						B.setElement(idx[1], idx[2], nIdx);
//						B.setElement(idx[2], idx[1], nIdx);
//					}
//					elements.add(new Triangle3(idx[0], idx[1], insertEleIdx[0], oldEle.label));
//					elements.add(new Triangle3(insertEleIdx[0], idx[2], idx[0], oldEle.label));
//				} else {
//					insertEleIdx[0] = eleIdx[2];
//					if (insertEleIdx[0] == 1) {
//						nIdx++;
//						insertX = 0.5 * (coord[0][0] + coord[2][0]);
//						insertY = 0.5 * (coord[0][1] + coord[2][1]);
//						nodes.add(new Node2D(insertX, insertY));
//						insertEleIdx[0] = nIdx;
//						B.setElement(idx[0], idx[2], nIdx);
//						B.setElement(idx[2], idx[0], nIdx);
//					}
//					elements.add(new Triangle3(idx[0], idx[1], insertEleIdx[0], oldEle.label));
//					elements.add(new Triangle3(insertEleIdx[0], idx[1], idx[2], oldEle.label));
//				}
//			} else {
//				elements.add(new Triangle3(idx, oldEle.label));
//			}
//		}
//
//		boundarys.clear();
//		int bCount = 0;
//		for (int i = 0; i < oldnb; i++) {
//			int[] idx = oldBoundarys.get(i).idx;
//			int val = B.getElement(idx[0], idx[1]);
//			if (val != 0) {
//				bCount++;
//			}
//		}
//
//		int newNb = oldnb + bCount;
//		nb = newNb;
//
//		double[][] coordB = new double[2][2];
//		for (int i = 0; i < oldnb; i++) {
//			Element ele = oldBoundarys.get(i);
//			for (int j = 0; j < 2; j++) {
//				coordB[j] = nodes.get(ele.idx[j]).coord.clone();
//			}
//			int val = B.getElement(ele.idx[0], ele.idx[1]);
//			if (val != 0) {
//				double[] newXY = bdFunc.action(coordB, ele.label, param);
//				nodes.set(val, new Node2D(newXY[0], newXY[1]));
//				boundarys.add(new Line2(ele.idx[0], val, ele.label));
//				boundarys.add(new Line2(val, ele.idx[1], ele.label));
//			} else {
//				boundarys.add(new Line2((Line2) (ele)));
//			}
//		}
//	}
//
//	@Override
//	public void setBoundarys(HashSet<UnorderIntegerArray> set) {
//		this.nb = set.size();
//		for(UnorderIntegerArray ele : set) {
//			this.boundarys.add(new Line2(ele.getIntSet(), 0));
//		}
//	}

	@Override
	public int getTecplotNt() {
		// TODO Auto-generated method stub
		return nt;
	}


	public int[][] getEdgesFromElement(int i) {
		int[][] result = new int[3][];
		int[] ele = elements[i];
		result[0] = new int[] {ele[0], ele[1]};
		result[1] = new int[] {ele[1], ele[2]};
		result[2] = new int[] {ele[2], ele[0]};
		return result;
	}

	public void toTecplotP0(String fileName, double[] sol) {
		assert (sol.length == nt);
		FileWriter fw = null;
		try {
			fw = new FileWriter(fileName);
			fw.write(Tecplot.solutionFileString + "\n");
			fw.write(Tecplot.dimVariableString[nDim - 1] + ", u\n");
			fw.write("ZONE NODES=" + 3 * nt + ", ELEMENTS=" + nt + ", DATAPACKING=POINT, ZONETYPE=" + getTecplotType());
			double[][] coord = new double[nPerElement][nDim];
			for (int i = 0; i < nt; i++) {
				int[] ele = elements[i];
				getCoordPerElement(ele, coord);
				for(int j = 0; j < nPerElement; j++) {
					for(int k = 0; k < nDim; k++) {
						fw.write(coord[j][k] + "  ");
					}
					fw.write(sol[i] + "\n");
				}
			}

			int idx = 0;
			for (int i = 0; i < nt; i++) {
				for (int j = 0, n = elements[i].length; j < n; j++) {
					fw.write((++idx) + "\t");
				}
				fw.write("\n");
			}
			fw.flush();
			fw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
}
